******************************************************
*******Engel curve parameter estimates ****************
******************************************************


*Set root data directory
local rootdir
cd "`rootdir'"



***Version 1 -- single unique base variety

local list grains veg
local agg regioncode round
local agg2 regioncode
local baseround=38
local otherround=66
local base x


use round`baseround'_edit, clear
append using round`otherround'_edit


**Appendix Figure A.7, group-level expenditure shares
preserve
cap: gen lnexp=log(totexp)
cap: gen sharegrains=xgrains/totexp
cap: gen shareveg=xveg/totexp
summ lnexp sharegrains shareveg if round==38, detail
set scheme s1mono
twoway (lowess sharegrains lnexp if lnexp>4.2 & lnexp<7.8 & sharegrains>0 & sharegrains<0.7 & round==38, xtitle(log expenditure) ytitle(Expenditure share)) (lowess shareveg lnexp if lnexp>4.2 & lnexp<7.8 & shareveg>0 & shareveg<0.17 & round==38)
graph export exp_shares.pdf, replace
restore


***lots to drop***
gen lnexp2=log(totexp-xfood)
drop *cloth* *bev* *intox* *light* *sugar* *proc* *meat* *oil* *milk* *food*
drop h*grains_* h*pulses_* h*veg_* h*fru_*
drop freemeal*
drop z mkt*
drop if totexp==0 | totexp==.


merge 1:1 round sector subround region subsample segment substratum fsu hhno using round`baseround'_hh, keepusing(mpce hhsize religion dmratioa dfratioa dmratios dfratios nic1 htype hgroup headed servant)
drop _merge
merge 1:1 round sector subround region subsample segment substratum fsu hhno using round`otherround'_hh, keepusing(mpce hhsize religion dmratioa dfratioa dmratios dfratios nic1 htype hgroup headed servant) update
drop _merge

merge m:m region round using mult`baseround', keepusing(regioncode)
drop _merge
merge m:m region using mult`otherround', keepusing(regioncode) update
drop _merge

drop region


drop if mpce>10000 | mpce<10
drop if hhsize>14 | hhsize==. | hhsize==0


**prices and price indexes (already done at aggregate levels?) - just merge in**
merge m:1 `agg' using pindex_r_revised
keep if _merge==3
drop _merge







***Select the base variety (and pick out its price too, then can drop quantities)

sort `agg'

foreach i in `list'{
gen base`i'=0
gen base`i'_alt=0
gen xbase`i'=0
gen xbase`i'_alt=0
gen aggbase`i'=0
gen basep`i'=0
levelsof n`i'tot, local(nmax)
egen hhbase`i'=rowmax(x`i'_1-x`i'_`nmax')

forvalues j=1(1)`nmax'{
by `agg2': egen totx=total(x`i'_`j'/qpindex`i')
by `agg': egen totx_alt=total(x`i'_`j'/qpindex`i')
replace base`i'=`j' if totx>xbase`i' & totx~=.
replace xbase`i'=totx if base`i'==`j'

replace base`i'_alt=`j' if totx_alt>xbase`i'_alt & totx_alt~=.
replace xbase`i'_alt=totx if base`i'==`j'

replace aggbase`i'=x`i'_`j' if base`i'==`j'
by `agg': egen p`i'=median((x`i'_`j'/q`i'_`j')*(qpindex`i'/pindex`i'))
replace basep`i'=p`i' if base`i'==`j'
drop totx p`i'
drop totx_alt
}
drop xbase`i' 
drop xbase`i'_alt

gen inbase`i'=0
replace inbase`i'=1 if aggbase`i'==hhbase`i' & aggbase`i'~=0
gen binbase`i'=0
replace binbase`i'=1 if aggbase`i'>0
gen cinbase`i'=0
replace cinbase`i'=1 if n`i'>0 & x`i'>0

bysort `agg2': egen minbase`i'=min(base`i'_alt)
bysort `agg2': egen maxbase`i'=max(base`i'_alt)
gen same`i'=0
replace same`i'=1 if minbase`i'==maxbase`i'
}
drop qgrains* qpulse* qveg* qfru*

drop if basegrains==0
summ base* *inbase*

summ same*
bysort regioncode: gen rcounter=_n
summ same* if rcounter==1


**97% of households consume some grains or veg
**92% of households consume the aggregate base variety
**78% (grains) and 52% (veg) consume the agg. base variety as their top variety


gen lnexp=log(mpce*hhsize)

bysort `agg': egen medmpce=median(mpce)
gen abovemed=0 if mpce~=.
replace abovemed=1 if mpce>medmpce & mpce~=.

***household demographic controls
gen lhhsize=log(hhsize)
gen rm=(dmratioa+dmratios)/hhsize
gen rf=(dfratioa+dfratios)/hhsize

***agriculture dummy
gen agg=0
replace agg=1 if nic1=="0"
replace agg=1 if sector==1 & (htype==1 | htype==4)

*hgroup
gen scst=1
replace scst=0 if hgroup==4 | hgroup==9
*religion
gen nonhindu=0
replace nonhindu=1 if religion==1

**pds
replace pds_qrice=0 if pds_qrice==.
replace pds_qwheat=0 if pds_qwheat==.

gen pdsqgrains=pds_qrice+pds_qwheat


cap: gen pds_xrice=0
cap: gen pds_xwheat=0
replace pds_xrice=0 if pds_xrice==.
replace pds_xwheat=0 if pds_xwheat==.
gen pdsxgrains=pds_xrice+pds_xwheat
gen pdssharegrains=pdsxgrains/xgrains

cap: gen pdsanygrains=0
cap: replace pdsanygrains=1 if pdsxgrains>0


gen pdsveg=0
gen pdsfru=0
gen pdspulses=0

**home share?
foreach i in `list'{
gen hshare`i'=hx`i'/x`i'
}


***village fixed effects...
egen cluster=group(round sector subround regioncode fsu)
bysort cluster: gen count=_N
drop if count<4 | count>10


****For eps,F regression
foreach i in `list'{
gen logx`i'=log(x`i')
gen logn`i'=log(n`i')
}


gen sigmagrains=2.16
gen sigmaveg=1.99


foreach i in `list'{
gen logpn`i'_agg=log((x`i'/aggbase`i')^(1/(sigma`i'-1)))
gen logpn`i'_hh=log((x`i'/hhbase`i')^(1/(sigma`i'-1)))
}




egen z=group(`agg')
egen maxz=max(z)
local maxz=maxz[1]


xtset cluster

gen year=1983
replace year=2009 if round==66

local pntype agg
local samp inbase
local clus cluster cluster
local iv i.headed 







****Appendix: Figure A.8********
preserve

local i grains
cap: gen share`i'=aggbase`i'/x`i'
cap: replace share`i'=1 if n`i'==1
cap: replace share`i'=log(share`i')
twoway (scatter share`i' logn`i' if `samp'`i'==1 & regioncode==141, xtitle("Log number of grain varieties") ytitle("Log exp. share of base grain") by(year, legend(off)) ) (lowess share`i' logn`i' if `samp'`i'==1 & regioncode==141)
graph save grains_141, replace

local i veg
cap: gen share`i'=aggbase`i'/x`i'
cap: replace share`i'=1 if n`i'==1
cap: replace share`i'=log(share`i')
twoway (scatter share`i' logn`i' if `samp'`i'==1 & regioncode==141, xtitle("Log number of veg. varieties") ytitle("Log exp. share of base veg.") by(year, legend(off)) ) (lowess share`i' logn`i' if `samp'`i'==1 & regioncode==141)
graph save veg_141, replace

graph combine grains_141.gph veg_141.gph, rows(2)
graph export psi_slope.png, replace


restore



summ lhhsize rm rf agg scst nonhindu pds* hshare* logn* logx* lnexp*

foreach i in `list'{
local controls lhhsize rm rf agg scst nonhindu pds*`i' hshare`i'

*a - ols
*c - IV

gen ec`i'_a=.
gen ec`i'_c=.

gen int`i'_a=.
gen int`i'_c=.

gen psi`i'_a=.
gen psi`i'_c=.

gen endog`i'_ec=.
gen endog`i'_psi=.

gen int`i'=logx`i'*abovemed

gen ecslopesplit`i'=.
gen ecslopesplit`i'_se=.

gen ecintsplit`i'=.
gen ecintsplit`i'_se=.

gen int`i'2=logn`i'*abovemed

gen psislopesplit`i'=.
gen psislopesplit`i'_se=.

forvalues j=1(1)`maxz'{
**OLS
capture: reg logn`i' logx`i' `controls' if z==`j' & `samp'`i'==1, vce(`clus')
capture: replace ec`i'_a=_b[logx`i'] if z==`j'
capture: replace int`i'_a=_b[_cons] if z==`j'


capture: reg logn`i' logx`i' abovemed int`i' `controls' if z==`j' & `samp'`i'==1, vce(`clus')
capture: replace ecslopesplit`i'=_b[int`i'] if z==`j'
capture: replace ecslopesplit`i'_se=_se[int`i'] if z==`j'
capture: replace ecintsplit`i'=_b[abovemed] if z==`j'
capture: replace ecintsplit`i'_se=_se[abovemed] if z==`j'


**IV
capture: ivregress 2sls logn`i' `controls' (logx`i'=`iv') if z==`j' & `samp'`i'==1, vce(`clus')
cap: estat endogenous 
cap: replace endog`i'_ec= r(p_regF) if z==`j' 
capture: replace ec`i'_c=_b[logx`i'] if z==`j'
capture: replace int`i'_c=_b[_cons] if z==`j'


capture: reg logpn`i'_`pntype' logn`i' `controls' if z==`j' & `samp'`i'==1, vce(`clus')
capture: replace psi`i'_a=_b[logn`i'] if z==`j'

capture: reg logpn`i'_`pntype' logn`i' abovemed int`i'2 `controls' if z==`j' & `samp'`i'==1, vce(`clus')
capture: replace psislopesplit`i'=_b[int`i'2] if z==`j'
capture: replace psislopesplit`i'_se=_se[int`i'2] if z==`j'

capture: ivregress 2sls logpn`i'_`pntype' `controls' (logn`i'=`iv') if z==`j' & `samp'`i'==1, vce(`clus')
cap: estat endogenous 
cap: replace endog`i'_psi= r(p_regF) if z==`j' 
capture: replace psi`i'_c=_b[logn`i'] if z==`j'

}

}


foreach i in `list'{
replace x`i'=. if x`i'==0
replace n`i'=. if n`i'==0
bysort `agg': egen x`i'_p50=median(x`i')
bysort `agg': egen x`i'_p90=pctile(x`i'), p(90)
bysort `agg': egen x`i'_p10=pctile(x`i'), p(10)

bysort `agg': egen n`i'_p50=median(n`i')
bysort `agg': egen n`i'_mean=mean(n`i')
bysort `agg': egen n`i'_p90=pctile(n`i'), p(90)
bysort `agg': egen n`i'_p10=pctile(n`i'), p(10)
}




bysort `agg' cluster: gen numhh=_N

foreach i in `list'{
levelsof n`i'tot, local(nmax)
gen rpresence`i'=0
gen vpresence`i'=0 if numhh>6
forvalues j=1(1)`nmax'{
bysort `agg': egen aggx`i'_`j'=total(x`i'_`j')
replace rpresence`i'=rpresence`i'+1 if aggx`i'_`j'>0 & aggx`i'_`j'~=.
drop aggx`i'_`j'
bysort `agg' cluster: egen aggx`i'_`j'=total(x`i'_`j')
replace vpresence`i'=vpresence`i'+1 if aggx`i'_`j'>0 & aggx`i'_`j'~=.
drop aggx`i'_`j'
}
bysort `agg': egen mean_vpresence`i'=mean(vpresence`i')
bysort `agg': egen median_vpresence`i'=median(vpresence`i')
}



duplicates drop `agg', force


foreach i in grains veg{
capture: gen ecsigslope`i'=abs(ecslopesplit`i'/ecslopesplit`i'_se)
capture: gen ecsigint`i'=abs(ecintsplit`i'/ecintsplit`i'_se)
capture: gen psisigslope`i'=abs(psislopesplit`i'/psislopesplit`i'_se)
}

 summ ecsigslope* ecsigint* psisigslope*, detail
bysort round: summ ecsigslope* ecsigint* psisigslope*, detail
bysort round: summ ecsigslope* ecsigint* psisigslope* if regioncode==141 
bysort round: summ ecsigslope* ecsigint* psisigslope* if regioncode==171

foreach i in grains veg{
foreach k in ecsigslope ecsigint psisigslope{
cap: gen sig5_`k'`i'=0 if `k'`i'~=.
cap: replace sig5_`k'`i'=1 if `k'`i'>1.96 & `k'`i'~=.
}

cap: gen sig5_ivec`i'=0 if endog`i'_ec~=.
cap: replace sig5_ivec`i'=1 if endog`i'_ec<=0.05

cap: gen sig5_ivpsi`i'=0 if endog`i'_psi~=.
cap: replace sig5_ivpsi`i'=1 if endog`i'_psi<=0.05

}



summ sig5_* endog*
bysort round: summ sig5_* endog*
bysort round: summ endog* if regioncode==141 | regioncode==171

local restrict if round==38

***Appendix Figure A.10 -- test for equality at mean

local i grains
twoway (kdensity ecsigslope`i' `restrict', title(`i') ytitle("Density of regions") xtitle("t-stat for null of parameter equality") xline(1.96) xlabel(1.96) legend(order(1 "EC slope" 2 "EC intercept" 3 "{&psi} slope"))) (kdensity ecsigint`i' `restrict') (kdensity psisigslope`i' `restrict')
graph save `i'_equality, replace
local i veg
twoway (kdensity ecsigslope`i' `restrict', title(`i') ytitle("Density of regions") xtitle("t-stat for null of parameter equality") xline(1.96) xlabel(1.96) legend(order(1 "EC slope" 2 "EC intercept" 3 "{&psi} slope"))) (kdensity ecsigint`i' `restrict') (kdensity psisigslope`i' `restrict')
graph save `i'_equality, replace
graph combine grains_equality.gph veg_equality.gph,  xcommon ycommon
graph export equality_rev1_1983.png, replace


****Appendix Figure A.11 -- Hausman specification test for IV estimates
*Hausman test
local i grains
twoway (kdensity endog`i'_ec `restrict', title(`i') ytitle("Density of regions") xtitle("p-value for Hausman test") xline(0.05) xlabel(0.05) legend(order(1 "EC equation" 2 "{&psi} equation"))) (kdensity endog`i'_psi `restrict')
graph save `i'_hausman, replace
local i veg
twoway (kdensity endog`i'_ec `restrict', title(`i') ytitle("Density of regions") xtitle("p-value for Hausman test") xline(0.05) xlabel(0.05) legend(order(1 "EC equation" 2 "{&psi} equation"))) (kdensity endog`i'_psi `restrict')
graph save `i'_hausman, replace
graph combine grains_hausman.gph veg_hausman.gph, xcommon ycommon
graph export hausman_rev1_1983.png, replace



*what else to keep? prices and such for decomposition? median expenditures and expenditure percentiles?
keep `agg' ec* int* psi* qpindex* basep* x*_p* n*_p* n*_mean endog*  *presence* base*

save parameters_rev1, replace






***Version 2 -- base variety varies by household [use aggregate price index]


local list grains veg 
local agg regioncode round
local agg2 regioncode
local baseround=38
local otherround=66
local base x




use round`baseround'_edit, clear
append using round`otherround'_edit

***lots to drop***
gen lnexp2=log(totexp-xfood)
drop *cloth* *bev* *intox* *light* *sugar* *proc* *meat* *oil* *milk* *food*
drop h*grains_* h*pulses_* h*veg_* h*fru_*
drop freemeal*
drop z mkt*
drop if totexp==0 | totexp==.


merge 1:1 round sector subround region subsample segment substratum fsu hhno using round`baseround'_hh, keepusing(mpce hhsize religion dmratioa dfratioa dmratios dfratios nic1 htype hgroup headed)
drop _merge
merge 1:1 round sector subround region subsample segment substratum fsu hhno using round`otherround'_hh, keepusing(mpce hhsize religion dmratioa dfratioa dmratios dfratios nic1 htype hgroup headed) update
drop _merge

merge m:m region round using mult`baseround', keepusing(regioncode)
drop _merge
merge m:m region using mult`otherround', keepusing(regioncode) update
drop _merge

drop region


drop if mpce>10000 | mpce<10
drop if hhsize>14 | hhsize==. | hhsize==0

**prices and price indexes (already done at aggregate levels?) - just merge in**
merge m:1 `agg' using pindex_r_revised
keep if _merge==3
drop _merge


drop *fru*


***Select the base variety (and pick out its price too, then can drop quantities)

sort `agg'

foreach i in `list'{
levelsof n`i'tot, local(nmax)
**pick each households base variety?
egen hhbase`i'=rowmax(x`i'_1-x`i'_`nmax')
gen basep`i'=qpindex`i'
gen cinbase`i'=0
replace cinbase`i'=1 if n`i'>0 & x`i'>0
}
drop qgrains* qpulse* qveg*

summ *inbase*, detail



gen lnexp=log(mpce*hhsize)




bysort `agg': egen medmpce=median(mpce)
gen abovemed=0 if mpce~=.
replace abovemed=1 if mpce>medmpce & mpce~=.




***household demographic controls
gen lhhsize=log(hhsize)
gen rm=(dmratioa+dmratios)/hhsize
gen rf=(dfratioa+dfratios)/hhsize

***agriculture dummy
gen agg=0
replace agg=1 if nic1=="0"
replace agg=1 if sector==1 & (htype==1 | htype==4)

*hgroup
gen scst=1
replace scst=0 if hgroup==4 | hgroup==9
*religion
gen nonhindu=0
replace nonhindu=1 if religion==1

**pds
replace pds_qrice=0 if pds_qrice==.
replace pds_qwheat=0 if pds_qwheat==.

gen pdsqgrains=pds_qrice+pds_qwheat

cap: gen pds_xrice=0
cap: gen pds_xwheat=0
replace pds_xrice=0 if pds_xrice==.
replace pds_xwheat=0 if pds_xwheat==.
gen pdsxgrains=pds_xrice+pds_xwheat
gen pdssharegrains=pdsxgrains/xgrains

cap: gen pdsanygrains=0
cap: replace pdsanygrains=1 if pdsxgrains>0

gen pdsveg=0
gen pdsfru=0
gen pdspulses=0

**home share
foreach i in `list'{
gen hshare`i'=hx`i'/x`i'
}


***village fixed effects...
egen cluster=group(round sector subround regioncode fsu)
bysort cluster: gen count=_N
drop if count<4 | count>10




****For eps,F regression
foreach i in `list'{
gen logx`i'=log(x`i')
gen logn`i'=log(n`i')
}



gen sigmagrains=2.16
gen sigmaveg=1.99


foreach i in `list'{
gen logpn`i'_hh=log((x`i'/hhbase`i')^(1/(sigma`i'-1)))
}




egen z=group(`agg')
egen maxz=max(z)
local maxz=maxz[1]

xtset cluster

local pntype hh
local samp cinbase
local clus cluster cluster
local iv i.headed 



summ lhhsize rm rf agg scst nonhindu pds* hshare* logn* logx* lnexp*

foreach i in `list'{
local controls lhhsize rm rf agg scst nonhindu pds*`i' hshare`i'

*a - ols
*c - IV

gen ec`i'_a=.
gen ec`i'_c=.

gen int`i'_a=.
gen int`i'_c=.

gen psi`i'_a=.
gen psi`i'_c=.

gen endog`i'_ec=.
gen endog`i'_psi=.

gen int`i'=logx`i'*abovemed

gen ecslopesplit`i'=.
gen ecslopesplit`i'_se=.

gen ecintsplit`i'=.
gen ecintsplit`i'_se=.

gen int`i'2=logn`i'*abovemed

gen psislopesplit`i'=.
gen psislopesplit`i'_se=.

forvalues j=1(1)`maxz'{

*OLS
capture: reg logn`i' logx`i' `controls' if z==`j' & `samp'`i'==1, vce(`clus')
capture: replace ec`i'_a=_b[logx`i'] if z==`j'
capture: replace int`i'_a=_b[_cons] if z==`j'


capture: reg logn`i' logx`i' abovemed int`i' `controls' if z==`j' & `samp'`i'==1, vce(`clus')
capture: replace ecslopesplit`i'=_b[int`i'] if z==`j'
capture: replace ecslopesplit`i'_se=_se[int`i'] if z==`j'
capture: replace ecintsplit`i'=_b[abovemed] if z==`j'
capture: replace ecintsplit`i'_se=_se[abovemed] if z==`j'

*IV
capture: ivregress 2sls logn`i' `controls' (logx`i'=`iv') if z==`j' & `samp'`i'==1, vce(`clus')
cap: estat endogenous 
cap: replace endog`i'_ec= r(p_regF) if z==`j' 
capture: replace ec`i'_c=_b[logx`i'] if z==`j'
capture: replace int`i'_c=_b[_cons] if z==`j'

capture: reg logpn`i'_`pntype' logn`i' `controls' if z==`j' & `samp'`i'==1, vce(`clus')
capture: replace psi`i'_a=_b[logn`i'] if z==`j'


capture: reg logpn`i'_`pntype' logn`i' abovemed int`i'2 `controls' if z==`j' & `samp'`i'==1, vce(`clus')
capture: replace psislopesplit`i'=_b[int`i'2] if z==`j'
capture: replace psislopesplit`i'_se=_se[int`i'2] if z==`j'

capture: ivregress 2sls logpn`i'_`pntype' `controls' (logn`i'=`iv') if z==`j' & `samp'`i'==1, vce(`clus')
cap: estat endogenous 
cap: replace endog`i'_psi= r(p_regF) if z==`j' 
capture: replace psi`i'_c=_b[logn`i'] if z==`j'

}

}


foreach i in `list'{
replace x`i'=. if x`i'==0
replace n`i'=. if n`i'==0
bysort `agg': egen x`i'_p50=median(x`i')
bysort `agg': egen x`i'_p90=pctile(x`i'), p(90)
bysort `agg': egen x`i'_p10=pctile(x`i'), p(10)

bysort `agg': egen n`i'_p50=median(n`i')
bysort `agg': egen n`i'_mean=mean(n`i')
bysort `agg': egen n`i'_p90=pctile(n`i'), p(90)
bysort `agg': egen n`i'_p10=pctile(n`i'), p(10)
}


bysort `agg' cluster: gen numhh=_N

foreach i in `list'{
levelsof n`i'tot, local(nmax)
gen rpresence`i'=0
gen vpresence`i'=0 if numhh>6
forvalues j=1(1)`nmax'{
bysort `agg': egen aggx`i'_`j'=total(x`i'_`j')
replace rpresence`i'=rpresence`i'+1 if aggx`i'_`j'>0 & aggx`i'_`j'~=.
drop aggx`i'_`j'
bysort `agg' cluster: egen aggx`i'_`j'=total(x`i'_`j')
replace vpresence`i'=vpresence`i'+1 if aggx`i'_`j'>0 & aggx`i'_`j'~=.
drop aggx`i'_`j'
}
bysort `agg': egen mean_vpresence`i'=mean(vpresence`i')
bysort `agg': egen median_vpresence`i'=median(vpresence`i')
}



duplicates drop `agg', force


foreach i in grains veg{
capture: gen ecsigslope`i'=abs(ecslopesplit`i'/ecslopesplit`i'_se)
capture: gen ecsigint`i'=abs(ecintsplit`i'/ecintsplit`i'_se)
capture: gen psisigslope`i'=abs(psislopesplit`i'/psislopesplit`i'_se)
}

 summ ecsigslope* ecsigint* psisigslope*, detail
bysort round: summ ecsigslope* ecsigint* psisigslope*, detail
bysort round: summ ecsigslope* ecsigint* psisigslope* if regioncode==141 
bysort round: summ ecsigslope* ecsigint* psisigslope* if regioncode==171

foreach i in grains veg{
foreach k in ecsigslope ecsigint psisigslope{
cap: gen sig5_`k'`i'=0 if `k'`i'~=.
cap: replace sig5_`k'`i'=1 if `k'`i'>1.96 & `k'`i'~=.
}

cap: gen sig5_ivec`i'=0 if endog`i'_ec~=.
cap: replace sig5_ivec`i'=1 if endog`i'_ec<=0.05

cap: gen sig5_ivpsi`i'=0 if endog`i'_psi~=.
cap: replace sig5_ivpsi`i'=1 if endog`i'_psi<=0.05

}


*what else to keep? prices and such for decomposition? median expenditures and expenditure percentiles?
keep `agg' ec* int* psi* qpindex* basep* x*_p* n*_p* n*_mean  endog* *presence*

save parameters_rev2, replace


summ ec* int* psi*, detail







*********Appendix Figure A.9 -- comparison of parameter estimates under Common Base (Version 1) or Household-specific base (Version 2)


use parameters_rev1, clear

local j a
foreach i in grains veg{
gen p`i'=basep`i'
gen eps`i'_`j'=psi`i'_`j'+(1/ec`i'_`j')
gen logF`i'_`j'=(log(psi`i'_`j'/eps`i'_`j'))-((eps`i'_`j'-psi`i'_`j')*int`i'_`j')-log(p`i')
gen F`i'_`j'=exp(logF`i'_`j')
}

foreach i of varlist epsgrains_a psigrains_a Fgrains_a epsveg_a psiveg_a Fveg_a ecgrains_a ecveg_a{
gen `i'_38=`i'
drop `i'
}
merge 1:1 regioncode round using parameters_rev2, update replace
foreach i in grains veg{
cap: drop p`i' logF`i'_`j'
gen p`i'=basep`i'
gen eps`i'_`j'=psi`i'_`j'+(1/ec`i'_`j')
gen logF`i'_`j'=(log(psi`i'_`j'/eps`i'_`j'))-((eps`i'_`j'-psi`i'_`j')*int`i'_`j')-log(p`i')
gen F`i'_`j'=exp(logF`i'_`j')
}


set scheme s1mono
twoway (scatter epsgrains_a epsgrains_a_38 if epsgrains_a>0 & epsgrains_a_38>0, title("Grains {&epsilon}") xtitle("Common base") ytitle("HH base") legend(off)) (lfit  epsgrains_a_38 epsgrains_a_38 if epsgrains_a_38>0)
graph save grains_eps, replace
twoway (scatter Fgrains_a Fgrains_a_38, title("Grains F") xtitle("Common base") ytitle("HH base") legend(off)) (lfit Fgrains_a_38 Fgrains_a_38)
graph save grains_F, replace
twoway (scatter psigrains_a psigrains_a_38, title("Grains {&psi}") xtitle("Common base") ytitle("HH base") legend(off)) (lfit psigrains_a_38 psigrains_a_38)
graph save grains_psi, replace

twoway (scatter epsveg_a epsveg_a_38 if epsveg_a>0 & epsveg_a_38>0, title("Veg {&epsilon}") xtitle("Common base") ytitle("HH base") legend(off)) (lfit epsveg_a_38 epsveg_a_38)
graph save veg_eps, replace
twoway (scatter Fveg_a Fveg_a_38, title("Veg F") xtitle("Common base") ytitle("HH base") legend(off)) (lfit Fveg_a_38 Fveg_a_38)
graph save veg_F, replace
twoway (scatter psiveg_a psiveg_a_38, title("Veg {&psi}") xtitle("Common base") ytitle("HH base") legend(off)) (lfit psiveg_a_38 psiveg_a_38)
graph save veg_psi, replace


graph combine grains_eps.gph grains_F.gph grains_psi.gph veg_eps.gph veg_F.gph veg_psi.gph, rows(2)
graph export parameters.pdf, replace










*****Feenstra CES on aggregate*********
use parameters_rev1, clear
duplicates drop regioncode, force
keep regioncode base*
save parameters_rev1_base, replace
************Aggregate Region*******

use round38_edit, clear
append using round66_edit

drop if totexp==0 | totexp==.
merge m:m region round using mult38, keepusing(regioncode)
drop _merge
merge m:m region round using mult66, keepusing(regioncode) update replace
duplicates drop

drop _merge
merge 1:1 round sector subround region subsample segment substratum fsu hhno using round38_hh, keepusing(mult) update replace
drop _merge
merge 1:1 round sector subround region subsample segment substratum fsu hhno using round66_hh, keepusing(mult) update replace
drop _merge

cap: gen mult=1
keep xveg* xgrains* nvegtot ngrainstot round sector subround regioncode subsample segment substratum fsu hhno mult

merge m:1 regioncode round using pindex_r_revised
drop _merge


merge m:1 regioncode using parameters_rev1_base
egen vill=group(round regioncode)


gen sigmagrains=2.16
gen sigmaveg=1.99

foreach i in grains veg{
levelsof n`i'tot, local(nmax)
bysort vill: egen aggx`i'=total((x`i'*mult)/pindex`i')
forvalues j=1(1)`nmax'{
by vill: egen aggx`i'_`j'=total((x`i'_`j'*mult)/pindex`i')
drop x`i'_`j'
}
}

duplicates drop vill, force 


preserve
keep if round==66


foreach i in grains veg{
levelsof n`i'tot, local(nmax)
ren aggx`i' paggx`i'
*ren nvill`i' pnvill`i'
forvalues j=1(1)`nmax'{
ren aggx`i'_`j' paggx`i'_`j'
}
}


save reg66, replace
restore

keep if round==38
joinby regioncode using reg66



foreach i in grains veg{
levelsof n`i'tot, local(nmax)
gen commonx`i'=0
gen pcommonx`i'=0

gen commonbasex`i'=0
gen pcommonbasex`i'=0

forvalues j=1(1)`nmax'{
gen common`i'_`j'=1 if aggx`i'_`j'>0 & paggx`i'_`j'>0
replace commonx`i'=commonx`i'+aggx`i'_`j' if common`i'_`j'==1
replace pcommonx`i'=pcommonx`i'+paggx`i'_`j' if common`i'_`j'==1

replace commonbasex`i'=aggx`i'_`j' if base`i'==`j'
replace pcommonbasex`i'=paggx`i'_`j' if base`i'==`j'

}
gen lambda`i'=1-((commonx`i'/aggx`i')/(pcommonx`i'/paggx`i'))^(1/(sigma`i'-1))
gen baselambda`i'=1-((commonbasex`i'/aggx`i')/(pcommonbasex`i'/paggx`i'))^(1/(sigma`i'-1))
bysort regioncode: egen meanlambda`i'=mean(lambda`i')
bysort regioncode: egen medianlambda`i'=median(lambda`i')
bysort regioncode: egen meanbaselambda`i'=mean(baselambda`i')
bysort regioncode: egen medianbaselambda`i'=median(baselambda`i')

gen daggx`i'=aggx`i'-paggx`i'
gen lndaggx`i'=ln(aggx`i')-ln(paggx`i')
bysort regioncode: egen meandaggx`i'=mean(daggx`i')
bysort regioncode: egen meanlndaggx`i'=mean(lndaggx`i')
}

duplicates drop regioncode, force


keep regioncode *lambda* mean*
save cesreg_3866, replace

summ *lambda*



****Village****

use round38_edit, clear
egen extra=group(region subround sector subsample segment substratum)
append using round66_edit
replace extra=0 if extra==.
drop if totexp==0 | totexp==.
merge m:m region round using mult38, keepusing(regioncode)
drop _merge
merge m:m region round using mult66, keepusing(regioncode) update replace

cap: gen mult=1
keep xveg* xgrains* nvegtot ngrainstot round sector subround regioncode subsample segment substratum fsu hhno mult extra

merge m:1 regioncode round using pindex_r_revised
drop _merge


merge m:1 regioncode using parameters_rev1_base


duplicates drop


egen vill=group(round fsu extra)
bysort vill: gen hcounter=_N
bysort round: tab hcounter
drop if hcounter>10
drop if hcounter<7
drop hcounter


gen sigmagrains=2.16
gen sigmaveg=1.99

summ pindex* x*

foreach i in grains veg{
levelsof n`i'tot, local(nmax)
bysort vill: egen aggx`i'=total(x`i'/pindex`i')
forvalues j=1(1)`nmax'{
by vill: egen aggx`i'_`j'=total(x`i'_`j'/pindex`i')
drop x`i'_`j'
}
}

duplicates drop vill, force 


preserve
keep if round==66


foreach i in grains veg{
levelsof n`i'tot, local(nmax)
ren aggx`i' paggx`i'
forvalues j=1(1)`nmax'{
ren aggx`i'_`j' paggx`i'_`j'
}
}

keep regioncode pagg*

save vill66, replace
restore

keep if round==38

egen zz=group(regioncode)
summ zz

joinby regioncode using vill66



foreach i in grains veg{
levelsof n`i'tot, local(nmax)
gen commonx`i'=0
gen pcommonx`i'=0

gen commonbasex`i'=0
gen pcommonbasex`i'=0
forvalues j=1(1)`nmax'{
gen common`i'_`j'=1 if aggx`i'_`j'>0 & paggx`i'_`j'>0
replace commonx`i'=commonx`i'+aggx`i'_`j' if common`i'_`j'==1
replace pcommonx`i'=pcommonx`i'+paggx`i'_`j' if common`i'_`j'==1

replace commonbasex`i'=aggx`i'_`j' if base`i'==`j'
replace pcommonbasex`i'=paggx`i'_`j' if base`i'==`j'

}
gen lambda`i'=1-((commonx`i'/aggx`i')/(pcommonx`i'/paggx`i'))^(1/(sigma`i'-1))
gen baselambda`i'=1-((commonbasex`i'/aggx`i')/(pcommonbasex`i'/paggx`i'))^(1/(sigma`i'-1))
bysort regioncode: egen meanlambda`i'=mean(lambda`i')
bysort regioncode: egen medianlambda`i'=median(lambda`i')
bysort regioncode: egen meanbaselambda`i'=mean(baselambda`i')
bysort regioncode: egen medianbaselambda`i'=median(baselambda`i')

gen daggx`i'=aggx`i'-paggx`i'
gen lndaggx`i'=ln(aggx`i')-ln(paggx`i')
bysort regioncode: egen meandaggx`i'=mean(daggx`i')
bysort regioncode: egen meanlndaggx`i'=mean(lndaggx`i')
}





duplicates drop regioncode, force

summ mean* median* 
summ mean* median* if regioncode==141
summ mean* median* if regioncode==171


keep regioncode *lambda*
save cesvill_3866, replace

summ *lambda*
















********Combining to do decomposition and welfare for ECV model**********

local list grains veg



use parameters_rev2, clear
**Use this code to generate the common base estimates displayed in Table A.9
*use parameters_rev1, clear


drop if xgrains_p50==0
summ psi*

sort regioncode round

foreach i in `list'{
gen p`i'=basep`i'

foreach j in a c{


gen eps`i'_`j'=psi`i'_`j'+(1/ec`i'_`j')

gen logF`i'_`j'=(log(psi`i'_`j'/eps`i'_`j'))-((eps`i'_`j'-psi`i'_`j')*int`i'_`j')-log(p`i')
gen F`i'_`j'=exp(logF`i'_`j')

*for COL later
gen eta`i'_`j'=psi`i'_`j'/eps`i'_`j'
gen const`i'_`j'= ( (eta`i'_`j'^(eta`i'_`j'/(1-eta`i'_`j'))) - (eta`i'_`j'^(1/(1-eta`i'_`j'))) )^(eta`i'_`j'-1)

forvalues k=10(40)90{
gen u0`i'_`k'_`j'=(x`i'_p`k'/(p`i' * (F`i'_`j'^(eta`i'_`j')) * const`i'_`j') )^(1/(1-eta`i'_`j'))
}


*predicted n
gen n`i'_p50_`j'= ((x`i'_p50*psi`i'_`j')/(eps`i'_`j'*F`i'_`j'*p`i'))^(1/(eps`i'_`j'-psi`i'_`j'))

}
}



gen base=1 if round==38
foreach i in `list'{

by regioncode: egen bp`i'=min(p`i'*base)
by regioncode: egen bx`i'_p50=min(x`i'_p50*base)
by regioncode: egen bx`i'_p10=min(x`i'_p10*base)
by regioncode: egen bx`i'_p90=min(x`i'_p90*base)

by regioncode: egen bn`i'_p50=min(n`i'_p50*base)
by regioncode: egen bn`i'_p10=min(n`i'_p10*base)
by regioncode: egen bn`i'_p90=min(n`i'_p90*base)

by regioncode: egen bn`i'_mean=min(n`i'_mean)

gen relx`i'_p50=log(x`i'_p50/bx`i'_p50)-log(p`i'/bp`i')

foreach j in a c{

foreach k in psi`i' eps`i' F`i' n`i'_p50 u0`i'_10 u0`i'_50 u0`i'_90{
by regioncode: egen b`k'_`j'=min(`k'_`j'*base)
}

*this will give welfare holding psi constant (only due to change in epsilon)
gen beta`i'_`j'=bpsi`i'_`j'/eps`i'_`j'
gen bconst`i'_`j'=( (beta`i'_`j'^(beta`i'_`j'/(1-beta`i'_`j'))) - (beta`i'_`j'^(1/(1-beta`i'_`j'))) )^(beta`i'_`j'-1)

gen reln`i'_p50_`j'=log(n`i'_p50_`j'/bn`i'_p50_`j')
gen relx`i'_p50_`j'=relx`i'_p50*(1/(eps`i'_`j'-psi`i'_`j'))
gen relpsi`i'_p50_`j'=log(psi`i'_`j'/bpsi`i'_`j')
replace relpsi`i'_p50_`j'=relpsi`i'_p50_`j'*(1/(eps`i'_`j'-psi`i'_`j'))
gen logx0`i'_`j'=log((bx`i'_p50*bpsi`i'_`j')/(bF`i'_`j'*bp`i'*beps`i'_`j'))
replace relpsi`i'_p50_`j'=relpsi`i'_p50_`j'+(logx0`i'_`j'  * ( (1/(eps`i'_`j'-psi`i'_`j')) -  (1/(eps`i'_`j'-bpsi`i'_`j')) ) )
gen relvc`i'_p50_`j'=log(beps`i'_`j'/eps`i'_`j')+log(bF`i'_`j'/F`i'_`j')
replace relvc`i'_p50_`j'=relvc`i'_p50_`j'*(1/(eps`i'_`j'-psi`i'_`j'))
replace relvc`i'_p50_`j'=relvc`i'_p50_`j'+( logx0`i'_`j' * ( (1/(eps`i'_`j'-bpsi`i'_`j')) -  (1/(beps`i'_`j'-bpsi`i'_`j')) ) )


foreach k in 10 50 90{
gen newx`i'_`k'_`j'=(bu0`i'_`k'_`j'^(1-eta`i'_`j')) * p`i' * (F`i'_`j'^(eta`i'_`j')) * const`i'_`j'
gen ratio`i'_`k'_`j'=newx`i'_`k'_`j'/bx`i'_p`k'
gen vargain`i'_`k'_`j'=1-(ratio`i'_`k'_`j'*(bp`i'/p`i'))

gen bnewx`i'_`k'_`j'=(bu0`i'_`k'_`j'^(1-beta`i'_`j')) * p`i' * (F`i'_`j'^(beta`i'_`j')) * bconst`i'_`j'
gen bratio`i'_`k'_`j'=bnewx`i'_`k'_`j'/bx`i'_p`k'
gen bvargain`i'_`k'_`j'=1-(bratio`i'_`k'_`j'*(bp`i'/p`i'))

}





gen vargdiff`i'_`j'=vargain`i'_90_`j'-vargain`i'_10_`j'
gen bvargdiff`i'_`j'=bvargain`i'_90_`j'-bvargain`i'_10_`j'


gen ces`i'_`j'=1- ((n`i'_p50^(-psi`i'_`j'))/(bn`i'_p50^(-bpsi`i'_`j')))
gen bces`i'_`j'=1- ((n`i'_p50^(-bpsi`i'_`j'))/(bn`i'_p50^(-bpsi`i'_`j')))

gen mces`i'_`j'=1- ((n`i'_mean^(-psi`i'_`j'))/(bn`i'_mean^(-bpsi`i'_`j')))
gen mbces`i'_`j'=1- ((n`i'_mean^(-bpsi`i'_`j'))/(bn`i'_mean^(-bpsi`i'_`j')))


}



local k 50
gen U0`i'=bx`i'_p`k'/(bp`i'*(bn`i'_p`k'^-bpsi`i'_a)) - bF`i'_a*(bn`i'_p`k'^beps`i'_a)
gen xhypeq`i'=(U0`i'+ F`i'_a*(bn`i'_p`k'^eps`i'_a))*(bp`i'*(bn`i'_p`k'^-bpsi`i'_a))
gen U1`i'=bx`i'_p`k'/(bp`i'*(bn`i'_p`k'^-bpsi`i'_a)) - F`i'_a*(bn`i'_p`k'^eps`i'_a)
gen fracxhypeq`i'=xhypeq`i'/bx`i'_p`k'
gen denom`i'=(bp`i'*(bn`i'_p`k'^-bpsi`i'_a))

gen fracgains`i'=(1-fracxhypeq`i')/bvargain`i'_`k'_a

}

summ fracxhypeq* if round==66, detail
summ fracgains* if round==66, detail




preserve

keep if round==66
keep regioncode relx*_p50 vargain* bvargain*
save ecv_comp3866, replace

restore






preserve
keep if round==66
keep regioncode vargain*_50_a bvargain*_50_a reln*_p50_a relx*_p50_a relpsi*_p50_a relvc*_p50_a vargdiff*_a bvargdiff*_a *ces*_a
save ecv_results3866, replace
restore



**variance decomposition
foreach i in `list'{
foreach j in a c{

gen vall`i'_`j'=0
foreach k in x psi vc{
bysort round: egen vrel`k'`i'_p50_`j'=sd(rel`k'`i'_p50_`j')
replace vall`i'_`j'=vall`i'_`j'+(vrel`k'`i'_p50_`j'^2)
}
foreach k in x psi vc{
replace vrel`k'`i'_p50_`j'=(vrel`k'`i'_p50_`j'^2)/vall`i'_`j'
}

}
}





***Many possible counter-factuals: 
*fixed costs equal to zero [consume all varieties] --> likely overstating magnitude of fixed costs...
*fixed costs equal to zero [consume same varieties]


***Welfare if setting fixed costs to zero?
gen nmaxgrains=18
gen nmaxveg=29

foreach i in grains veg{
foreach j in 10 50 90{
*no fixed cost all n
gen refx`i'_`j'=(u0`i'_`j'_a^(1-eta`i'_a)) * p`i'* (F`i'_a^(eta`i'_a)) * const`i'_a

gen refn`i'_p`j'=((x`i'_p50*psi`i'_a)/(eps`i'_a*F`i'_a*p`i'))^(1/(eps`i'_a-psi`i'_a))

gen eqx`i'_`j'_noF=u0`i'_`j'_a* p`i'* (nmax`i'^(-psi`i'_a))

*no fixed cost fixed n
gen eqx`i'_`j'_noFn=u0`i'_`j'_a* p`i'* (refn`i'_p`j'^(-psi`i'_a))

*fixed cost that leads to consuming all varieties [at each percentile]
gen F_hyp`i'_`j'= ((x`i'_p`j'/p`i')*(psi`i'_a/eps`i'_a))/(nmax`i'^(eps`i'_a-psi`i'_a))
gen eqx`i'_`j'_F=(u0`i'_`j'_a^(1-eta`i'_a)) * p`i' * (F_hyp`i'_`j'^(eta`i'_a)) * const`i'_a

*fixed cost that leads to consuming all varieties [for 10th percentile -> extra infra-marginal gain for richer guys]
gen eqx`i'_`j'_Fmin=(u0`i'_`j'_a^(1-eta`i'_a)) * p`i' * (F_hyp`i'_10^(eta`i'_a)) * const`i'_a
}
}



foreach i in grains veg{
foreach z in eqx`i'_50_noF eqx`i'_50_noFn eqx`i'_50_F eqx`i'_50_Fmin{
replace `z'=1 - (`z'/x`i'_p50)
}
summ eqx`i'_50_noF eqx`i'_50_noFn eqx`i'_50_F eqx`i'_50_Fmin if round==66 &  eqx`i'_50_Fmin~=0, detail
}



keep if round==66
merge 1:1 regioncode using region_data
drop _merge
merge 1:1 regioncode using cesreg_3866, keepusing(lambdagrains lambdaveg)
replace lambdagrains=-lambdagrains
replace lambdaveg=-lambdaveg

drop _merge
merge 1:1 regioncode using cesvill_3866, keepusing(medianlambdagrains medianlambdaveg)
replace medianlambdagrains=-medianlambdagrains
replace medianlambdaveg=-medianlambdaveg
drop _merge


local regspec a
foreach k in grains veg{
pwcorr `j' reln`k'_p50_`regspec' relx`k'_p50_`regspec' relvc`k'_p50_`regspec' bvargain`k'_50_`regspec' vargain`k'_50_`regspec' bvargdiff`k'_`regspec' medianlambda`k' density43, sig
pwcorr `j' reln`k'_p50_`regspec' relx`k'_p50_`regspec' relvc`k'_p50_`regspec' bvargain`k'_50_`regspec' vargain`k'_50_`regspec' bvargdiff`k'_`regspec' medianlambda`k' lnpce43, sig
}

pwcorr density43 lnpce43, sig

local regspec a
local qvar qtiledensity43

foreach k in grains veg{
foreach j in reln`k'_p50_`regspec' relx`k'_p50_`regspec' relpsi`k'_p50_`regspec' relvc`k'_p50_`regspec' bvargain`k'_50_`regspec' vargain`k'_50_`regspec' bvargdiff`k'_`regspec' lambda`k' medianlambda`k'{

replace `j'=. if `j'==0
summ `j'
egen mean_`j'=mean(`j')
egen q1_`j'=mean(`j'*`qvar'_1)
egen q5_`j'=mean(`j'*`qvar'_5)
egen q3_`j'=mean(`j'*`qvar'_3)
}
}
keep mean_* q1_* q5_* q3_*
drop mean*presence*
duplicates drop

expand 9
gen meangrains=.
gen meanveg=.
gen q1grains=.
gen q1veg=.
gen q5grains=.
gen q5veg=.
gen q3grains=.
gen q3veg=.
gen row=""
replace row="Change in variety ($%\Delta n$)" if _n==1
replace row="Expenditure component (X/b)" if _n==2
replace row="Int. margin component ($\psi$)" if _n==3
replace row="Variety cost component ($ F, \epsilon$)" if _n==4
replace row="Change in variety cost ($ F, \epsilon$) only" if _n==5
replace row="Change in variety cost and int. margin ($ F, \epsilon, \psi$) " if _n==6
replace row="Rich bias: Gain at 90th vs. 10th percentile ($ F, \epsilon$)" if _n==7
*replace row="Rich bias: Gain at 90th vs. 10th percentile ($F, \epsilon, \psi$)" if _n==8
replace row="Region aggregate CES" if _n==8
replace row="Village aggregate CES (Region median)" if _n==9




foreach k in mean q1 q3 q5{
foreach i in grains veg{
ren `k'_bvargdiff`i'_a `k'_bvargdiff`i'_p50_a
ren `k'_bvargain`i'_50_a `k'_bvargain`i'_p50_a
ren `k'_vargain`i'_50_a `k'_vargain`i'_p50_a
}
}

foreach k in mean q1 q3 q5{
foreach i in grains veg{
local counter=1
foreach j in reln relx relpsi relvc bvargain vargain bvargdiff {
replace `k'`i'=`k'_`j'`i'_p50_a if _n==`counter'
local counter=`counter'+1
}
foreach j in lambda medianlambda{
replace `k'`i'=`k'_`j'`i' if _n==`counter'
local counter=`counter'+1
}
}
}

keep row meangrains q1grains q5grains meanveg q1veg q5veg
order row meangrains q1grains q5grains meanveg q1veg q5veg
browse



